Algorithm xxx — ORTHPOL: A package 

of routines for generating orthogonal 
polynomials and Gauss-type quadrature rules^ 
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ABSTRACT. A collection of subroutines and examples of their uses, as well as 
the underlying numerical methods, are described for generating orthogonal poly- 
nomials relative to arbitrary weight functions. The object of these routines is to 
produce the coefficients in the three-term recurrence relation satisfied by the or- 
thogonal polynomials. Once these are known, additional data can be generated, 
such as zeros of orthogonal polynomials and Gauss- type quadrature rules, for which 
routines are also provided. 
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1. INTRODUCTION 

Classical orthogonal polynomials, such as those of Legendre, Chebyshev, 
Laguerre and Hermite, have been used for purposes of approximation in 
widely different disciplines and over a long period of time. Their popularity 
is due in part to the ease with which they can be employed, and in part 
to the wealth of analytic results known for them. Widespread use of non- 
classical orthogonal polynomials, in contrast, has been impeded by a lack of 
effective and generally applicable constructive methods. The present set of 
computer routines has been developed over the past ten years in the hope 
of remedying this impediment and of encouraging the use of nonstandard 
orthogonal polynomials. A number of applications indeed have already been 
made, for example to numerical quadrature (Cauchy principal value integrals 
with coth-kernel [38], Hilbert transform of Jacobi weight functions [37], inte- 
gration over half-infinite intervals [28], rational Gauss-type quadrature [30, 
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31]), to moment-preserving spline approximation [21,35,11], to the summa- 
tion of slowly convergent series [26,27], and, perhaps most notably, to the 
proof of the Bieberbach conjecture [24]. 

In most applications, orthogonality is with respect to a positive weight 
function, w, on a given interval or union of intervals, or with respect to 
positive weights, Wi, concentrated on a discrete set of points, {xj}, or a 
combination of both. For convenience of notation, we subsume all these 
cases under the notion of a positive measure, dX, on the real line R. That 
is, the respective inner product is written as a Riemann-Stieltjes integral. 



where the function X{t) is the indefinite integral of w for the continuous 
part, and a step function with jumps Wi at for the discrete part. We 
assume that (1.1) is meaningful whenever u,v are polynomials. There is 
then defined a unique set of (monic) orthogonal polynomials, 



We speak of "continuous" orthogonal polynomials if the support of dX is 
an interval, or a union of intervals, of "discrete" orthogonal polynomials 
if the support of dX consists of a discrete set of points, and of orthogonal 
polynomials of "mixed type" if the support of dX has both a continuous and 
discrete part. In the first and last case, there are infinitely many orthogonal 
polynomials — one for each degree — , whereas in the second case there 
are exactly N orthogonal polynomials, ttq, tti, . . . , ttjv-i, where N is the 
number of support points. In all cases, we denote the polynomials by vr^.(-) = 
TTki • ',dX), or 7rfe( • ;u)), if we want to indicate their dependence on 
the measure dX or weight function w, and use similar notations for other 
quantities depending on dX or w. 

It is a distinctive feature of orthogonal polynomials, compared to other 
orthogonal systems, that they satisfy a three-term recurrence relation, 




(1.1) 



T^k{t) = 1^^+ lower-degree terms, k = 0,1,2, ... , 



(1.2) 



(7rfc,7r£) = if k=>^i. 



vrfe+i(t) = {t- ak)TTk{t) - PkT^k-iit), k = 0,1,2,... 



(1.3) 



7ro(t) = l, 7r_i(i) = 0, 
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with coefficients ak = ak{dX) G R, = Pk{dX) > that are uniquely 
determined by the measure dX. By convention, the coefficient /3o, which 
multiphes 7r_i = in (1.3), is defined by 



The knowledge of these coefficients is absolutely indispensable for any sound 
computational use and application of orthogonal polynomials [19,25]. One 
of the principal objectives of the present package is precisely to provide 
routines for generating these coefficients. Routines for related quantities are 
also provided, such as Gauss-type quadrature weights and nodes, hence also 
zeros of orthogonal polynomials. 

Occasionally (e.g., in [21,35,11,30,31]), one needs to deal with indefinite 
(i.e., sign-changing) measures dX. The positivity of the Pk is then no longer 
guaranteed, indeed not even the existence of all orthogonal polynomials. 
Nevertheless, our methods can still be formally applied, albeit at the risk of 
possible breakdowns or instabilities. 

There are basically four methods used here to generate recursion coef- 
ficients: (i) Methods based on explicit formulae. These relate to classical 
orthogonal polynomials, and arc implemented in the routine recur of §2. 
(ii) Methods based on moment information. These are dealt with in §3, and 
are represented by a single routine, cheb. Its origin can be traced back to 
work of Chebyshev on discrete least squares approximation, (iii) Bootstrap 
methods based on inner product formulae for the coefficients, and orthogonal 
reduction methods. We have attributed the idea for the former method to 
Stieltjes and referred to it in [19] as the Stieltjes procedure. The prototype 
is the routine sti in §4, applicable for discrete orthogonal polynomials. An 
alternative routine is lancz, which accomplishes the same purpose, but uses 
the method of Lanczos. Either of these routines can be used in mcdis, which 
applies to continuous as well as mixed-type orthogonal polynomials. In con- 
trast to all previous routines, mcdis uses a discretization process, and thus 
furnishes only approximate answers whose accuracies can be controlled by 
the user. The routine, however, is by far the most sophisticated and flexi- 
ble routine in this package, one that requires, or can greatly benefit from, 
ingenuity of the user. The same kind of discretization is also applicable to 
moment-related methods, yielding the routine mccheb. (iv) Modification al- 
gorithms. These are routines generating recursion coefficients for measures 
modified by a rational factor, utilizing the recursion coefficients of the orig- 
inal measure, which are assumed to be known. They can be thought of as 




(1.4) 



3 



algorithmic implementations of the Christoffel, or generalized Christoffel, 
theorem and are incorporated in the routines chri and gchri of §5. An 
important application of all these routines is made in §6, where routines are 
provided that generate the weights and nodes of quadrature rules of Gauss, 
Gauss-Radau, and Gauss-Lobatto type. 

Each routine has a single-precision and double-precision version with 
similar names, except for the prefix "d" in double-precision procedures. The 
latter are generally a straightforward translation of the former. An excep- 
tion is the routine dlga used in drecur for computing the logarithm of the 
gamma function, which employs a different method than the single-precision 
companion routine alga. 

All routines of the package have been checked for ANSI conformance, 
and have been tested on two computers: the Cyber 205 and a Sun 4/670 
MP workstation. The former has machine precisions 7.11 X 10~^^, 

e'^ 5.05 X 10~^^ in single and double precision, respectively, while the 
latter has ^ 5.96 x 10"^, e"* « 1.11 x 10"^^. The Cyber 205 has a 
large floating-point exponent range, extending from approximately —8617 
to -1-8645 in single as well as double precision, whereas the Sun 4/670 has 
the rather limited exponent range [-38, 38] in single precision, but a larger 
range, [-308, 308], in double precision. All output cited relates to work on 
the Cyber 205. 

The package is organized as follows: Section contains (slightly amended) 
netlib routines, namely rlmach and dlmach, providing basic machine con- 
stants for a variety of computers. Section 1 contains all the driver routines 
— named testl, test2, etc. — which are used (and described in the body 
of this paper) to test the subroutines of the package. The complete output of 
each test is listed immediately after the driver. Sections 2-6 constitute the 
core of the package: the single- and double-precision subroutines described 
in the equally-numbered sections of this paper. All single-precision routines 
are provided with comments and instructions for their use. These, of course, 
apply to the double-precision routines as well. 

2. CLASSICAL WEIGHT FUNCTIONS 

Among the most frequently used orthogonal polynomials are the Jacobi 
polynomials, generalized Laguerre polynomials, and Hermite polynomials, 
supported respectively on a finite interval, half-infinite interval, and the 



4 



whole real line. The respeetive weight funetions are 

=t(;("./3)(t) = + on (-1,1), a>-l,P>-l: Jacobi, 

(2.1) 

w{t) = w^'^\t) = f^e * on (0,00), a > — 1 : generalized Laguerre, (2.2) 
u'(t) = e^*^ on (—00,00): Hermite. (2.3) 

Special cases of the Jacobi polynomials are the Legendre polynomials {a = 
P = 0), the Chebyshev polynomials of the first (a = /? = —5), second 
(a = /3 = ^), third (a = — /3 = — ^) and fourth {a = —(3 = ^) kind, and the 
Gegenbauer polynomials (a = /3 = A— ^). The Laguerre polynomials are 
the special case a = of the generalized Laguerre polynomials. 

For each of these polynomials, the corresponding recursion coefficients 
= (^k{w), Pk = Pk{w) arc explicitly known (sec, e.g., [5, pp. 217-221]), 
and are generated in single precision by the routine recur. Its calling se- 
quence is 

recur (n , ipoly , al , be , a , b , ierr ) . 

On entry, 

n is the number of recursion coefficients desired; 

type integer 

ipoly is an integer identifying the polynomial as 
follows: 

1 = Legendre polynomial on (—1,1) 

2 = Legendre polynomial on (0,1) 

3 = Chebyshev polynomial of the first kind 

4 = Chebyshev polynomial of the second kind 

5 = Chebyshev polynomial of the third kind 

6 = Jacobi polynomial with parameters al , be 

7 = generalized Laguerre polynomial with pa- 

rameter al 

8 = Hermite polynomial 

al.be are the input parameters a,/3 for Jacobi and 
generalized Laguerre polynomials; type real; 
they are only used if ipoly = 6 or 7, and in 
the latter case only al is used. 
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On return, 



are real arrays of dimension n with a.{k), b(A;) 
containing the coefficients a^-i, /3fe-i, respec- 
tively, = 1, 2, . . . , n 
is an error flag, where 
ierr = on normal return, 
ierr = 1 if either al or be are out of range 

when ipoly = 6 or ipoly=7, 
ierr = 2 if there is potential overflow in the 
evaluation of Pq when ipoly = 6 or 
ipoly = 7; in this case Pq is set equal 
to the largest machine-representable 
number, 
ierr = 3 if n is out of range, 
ierr = 4 if ipoly is not one of the 
admissible integers. 

No provision has been made for Chebyshev polynomials of the fourth kind, 
since their recursion coefficients are obtained from those for the third-kind 
Chebyshev polynomials simply by changing the sign of the a's (and leaving 
the /3's unchanged). 

The corresponding double-precision routine is drecur; it has the same 
calling sequence, except for real data types now being double-precision. 

In the cases of Jacobi polynomials (ipoly = 6) and generalized Laguerre 
polynomials (ipoly = 7), the recursion coefficient /3o (and only this one) 
involves the gamma function T. Accordingly, a function routine, alga, is 
provided that computes the logarithm In T of the gamma function, and a 
separate routine, gamma, computing the gamma function by exponentiating 
its logarithm. Their calling sequences are 

function alga(x) 
fvmction gamma (x, ierr) , 

where ierr is an output variable set equal to 2 or depending on whether 
the gamma function does, or does not, overflow, respectively. The corre- 
sponding double-precision routines have the names dlga and dgamma. All 
these routines require machine-dependent constants for reasons explained 
below. 



a,b 
ierr 
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The routine alga is based on a rational approximation valid on the 



interval 



where 



1 3 
2' 2 



Outside this interval, the argument x is written as 



Xe + m, 



Xp. = < 



x — [xj + 1 if X — [xj < ^ , 
X — [_x\ otherwise 



is in the interval | 



and m > — 1 is an integer. If m = — 1 (i.e., < x < 

^), then In r(x) = In r(xe)— In x, while for m > 0, one computes In r(x) = 
In r(xe) + In p, where p = Xe(xe + !)••• (xg + m — 1). If m is so large, say 
ITT' > iTT-o, that the product p would overflow, then In p is computed (at a 
price!) as In p = In Xg + In (xg + l) + • • • + In (xg + m — 1). It is here, where a 
machine-dependent integer is required, namely mo = smallest integer m such 
that 1 • 3 • 5 • • • (2m + l)/2'" is greater than or equal to the largest machine- 
representablc number, R. By Stirling's formula, the integer mo is seen to be 
the smallest integer m satisfying ((m + l)/e)ln((m + l)/e) > (Ini? — ^In8)/e, 
hence equal to [e • t{{\nR — ^ln8)/e)J, where t{y) is the inverse function of 
y = tint. For our purposes, the low-accuracy approximation of t(y), given 
in [16, pp. 51-52], and implemented in the routine t, is adequate. 

The rational approximation chosen on ^, | is one due to W.J. Cody 
and K.E. Hillstrom, namely the one labeled n = 7 in Table II of [9]. It 
is designed to yield about 16 correct decimal digits (cf. [9, Table I]), but 
because of numerical cancellation furnishes only about 13-14 correct decimal 
digits. 

Since rational approximations for In T having sufficient accuracies for 
double-precision computation do not seem to be available in the literature, 
we used a different approach for the routine dlga, namely the asymptotic 
approximation (cf. [1, Eq. 6.1.42], where the constants are Bernoulli 
numbers) 



Inr(y) = (y- i)lny-2/ + iln(27r) 



+ E 



m=l 



^^^^ y-i^rn-1) ^ j^^y^ 



(2.4) 
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for values of y > large enough to have 



\Rniy)\ < \ 10-^ (2.5) 

where d is the number of decimal digits carried in double-precision arith- 
metic, another machine-dependent real number. If (2.5) holds for y > yg, 
and if X > yo) we compute In r{x) from the asymptotic expression (2.4) 
(where y = X and the remainder term is neglected). Otherwise, we let ko be 
the smallest positive integer k such that x + k > yg, and use 

In r{x) = In r(x + ko) - In {x{x + 1) ■ ■ ■ {x + ko - 1)), (2.6) 

where the first term on the right is computed from (2.4) (with y = x + ko). 
Since for y > 0, 

(2n -I- 2)(2n -|- 1) 



(cf. [1, Eq. 6.1.42]), the inequality (2.5) is satisfied if 

1 



> exp 



2n + 1 



din 10 + In ^'^2"+^' 



(2n + 2)(2n + l) 



(2.7) 



In our routine dlga, we have selected n = 9. For double-precision ac- 
curacy on the Cyber 205, we have d « 28.3, for which (2.7) then gives 
y > exp{.121188 ... d+ .053905 32.6. 

For single-precision calculation we selected the method of rational ap- 
proximation, rather than the asymptotic formula (2.4) and (2.6), since wc 
found that the former is generally more accurate, by a factor, on the average, 
of about 20 and as large as 300. Neither method yields full machine accu- 
racy. The former, as already mentioned, loses accuracy in the evaluation of 
the approximation. The latter suffers loss of accuracy because of cancella- 
tion occurring in (2.6), which typically amounts to a loss of 2-5 significant 
decimal digits in the gamma function itself. 

Since these errors affect only the coefficient /?o (and only if ipoly = 6 
or 7), they are of no consequence unless the output of the routine recur 
serves as input to another routine, such as gauss (cf. §6), which makes 
essential use of /3o- In this case, for maximum single-precision accuracy, it 
is recommended that /3o be first obtained in double precision by means of 
drecur with n = 1, and then converted to single precision. 
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3. MOMENT-RELATED METHODS 



It is a well-known fact that the first n recursion coefficients ak{dX), 
Pk{dX), k = 0,1, . . . ,n — 1 (cf. (1.3)), are uniquely determined by the first 
2n moments Hk of the measure dX, 

^jij^ = ^k{d\)= [ t^d\{t), A; = 0,l,2,...,2n-1. (3.1) 

Formulae are known, for example, which express the a's and /3's in terms of 
Hankel determinants in these moments. The trouble is that these formulae 
become increasingly sensitive to small errors as n becomes large. There is an 
inherent reason for this: the underlying (nonlinear) map Kn' R^" R-^" 
has a condition number, cond Kn, which grows exponentially with n (cf. 
[19, §3.2]). Any method that attempts to compute the desired coefficients 
from the moments in (3.1), therefore, is doomed to fail, unless n is quite 
small, or extended precision is being employed. That goes in particular for 
an otherwise elegant method due to Chebyshev (who developed it for the 
case of discrete measures d\) that generates the a's and /3's directly from 
the moments (3.1), bypassing determinants altogether (cf. [4], [19, §2. 3]). 

Variants of Chebyshcv's algorithm with more satisfactory stability prop- 
erties have been developed by Sack and Donovan [46] and Wheeler [50] (in- 
dependently of Chebyshev's work). The idea is to forego the moments (3.1) 
as input data, and instead depart from so-called modified moments. These 
are defined by replacing the power in (3.1) by an appropriate polynomial 
Pk{t) of degree k, 

Uk = Uk{dX) = I Pk{t)dX{t), k = 0,1,2,..., 2n-l. (3.2) 
./ R 

For example, pk could be one of the classical orthogonal polynomials. More 
generally, we shall assume that {pk} are monic polynomials satisfying a 
three-term recurrence relation similar to the one in (1.3), 



Pk+i{t) = {t- ak)pkit) - bkPk-iit), k = 0,l,2,... , 

(3.3) 

Po{t) = l, p-i{t) = 0, 

with coefficients a*; G R, 6jk > that are known. (In the special case = 0, 
6fe = 0, we are led back to powers and ordinary moments.) There now 
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exists an algorithm, called modified Chebyshev algorithm in [19, §2.4], which 
takes the 2n modified moments in (3.2) and the 2n—l coefficients {ofcl^^o^, 
{^fc}fc=o^ in (3.3), and from them generates the n desired coefficients ak{dX), 
j3k{d\), /c = 0, 1, . . . , n — 1. It generalizes Chebyshcv's algorithm, which can 
be recovered (if need be) by putting = = 0. The modified Chebyshev 
algorithm is embodied in the subroutine cheb, which has the calling sequence 



cheb (n , a , b , f nu , alpha , beta , s , ierr , sO , s 1 , s2) 
dimension a(*) ,b(*) ,fnu(*) ,alpha(n) ,beta(n) ,s(n) , 
sO(*) ,sl(*) ,s2(*) 



On entry, 

n is the number of recursion coefficients desired; 

type integer 

a,b are arrays of dimension 2xn-l holding the 

coefficients a{k) = ak-i, b(/c) = bk-i, k = 
l,2,...,2n- 1 

f nu is an array of dimension 2 x n holding the mod- 

ified moments f nu(A;) = ffc-i, k = 1,2, . . . ,2x 



On return. 



alpha, beta are arrays of dimension n containing the de- 
sired recursion coefficients alpha (k) = ak-i, 
beta {k) = Pk-i, A; = 1, 2, . . . , n 

s is an array of dimension n containing the num- 

bers s(A;) = Jp^ Trl_idX, k = 1,2, . . . , n 

ierr is an error flag, equal to on normal return, 

equal to 1 if | i^o I is less than the machine zero, 
equal to 2 if n is out of range, equal to — {k— 1) 
if s(A;), /c = 1, 2, . . . , n, is about to underflow, 
and equal to -|-(A:— 1), if it is about to overflow. 
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The arrays sO, si, s2 of dimension 2xn are needed for working space. 



There is again a map Kn'. R R' 



underlying the modified Cheby- 



shev algorithm, namely the map taking the 2n modified moments into the 
n pairs of recursion coefficients. The condition of the map Kn (actually of a 
somewhat different, but closely related map) has been studied in [19, §3.3] 
and [23] in the important case where the polynomials pk defining the mod- 
ified moments arc themselves orthogonal polynomials, Pki') = Pk{ • ',dl^), 
with respect to a measure (for example, one of the classical ones) for 
which the recursion coefficients a^, bj- are known. The upshot of the anal- 
ysis then is that the condition of is characterized by a certain positive 
polynomial gf„ ( ■ ; dX) of degree 4n — 2, depending only on the target measure 
d\, in the sense that 



Thus, the numerical stability of the modified Chebyshev algorithm is deter- 
mined by the magnitude of gn on the support of d/j,. 

The occurrence of underflow [overflow] in the computation of the a's 
and /3's, especially on computers with limited exponent range, can often be 
avoided by multiplying all modified moments by a sufficiently large [small] 
scaling factor before entering the routine. On exit, the coefficient Po (and 
only this one!) then has to be divided by the same scaling factor. (There may 
occur harmless underflow of auxiliary quantities in the routine cheb, which 
is difficult to avoid since some of these quantitites actually are expected to 
be zero.) 

Example 3.1. dX^{t) = [(1 - u;'^f){l - t'^)]-'^/'^dt on (-1,1), < co < 1. 

This example is of some historical interest, in that it has already been 
considered by Christoffel [8, Example 6]; sec also [44]. Computationally, the 
example is of interest as there are empirical reasons to believe that for the 
choice diJ,{t) = (1 — t^)^^^'^dt on (—1,1) — which appears rather natural 
— the modified Chebyshev algorithm is exceptionally stable, uniformly in 
n, in the sense that in (3.4) one has Qn < 1 on supp dfi for all n (cf. [22, 
Example 5.2]). With the above choice of d^u, the polynomials pk are clearly 
the Chebyshev polynomials of the first kind, pq = Tq, pk = 2~^^~^^Tk, k > 1, 
and the modified moments are given by 




(3.4) 




(3.5) 
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They are expressible in terms of the Fourier coefficients Cr(a;^) in 

oo 

(1 - iJ sin2 0)-V2 = ^^(^2) ^ 2 ^ cos 2re (3.6) 

by means of (cf. [19, Example 3.3]) 



r=l 



^2m-l = 



(3.7) 



> m = l,2,3, . 



The Fourier coefficients {Cr(w^)} in turn can be accurately computed as 
minimal solution of a certain three-term recurrence relation (see [19, pp. 
310-311]). 

The ordinary moments 



k — 1,2,3,..., 



(3.^ 



likewise can be expressed in terms of the Fourier coefficients Cr{oj^) by writ- 
ing t^"^ as a linear combination of Chebyshev polynomials Tq, . . . ^T2m 
(cf. Eq. (22) on p. 454 of [43]). The result is 



At2 



22m 
/^2m-l = 0, 



-|Nm„ rn 

r E(-i)''7^^c/m-.(c.'), 

r=0 



m = 1, 2, 3, 



(3.9) 



where 



(m) 

7o 



7, 



(m) _ 2m + 1 - r (m) 



7r-i, r = 1,2, . . . ,m - 1, 



(3.10) 
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(m) _ 1 ^(m) 



2m 



7m-l- 



The driver routine testl (in §1 of the package) generates for = 
.1(.2).9, .99, .999 the first n recurrence coefficients Pk{d^u)) (all = 0), both 
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in single and double precision, using modified moments if niodinoin= . true . , 
and ordinary moments otherwise. In the former case, n = 80, in the latter, 
n = 20. It prints the double-precision values of (3^ together with the rela- 
tive errors of the single-precision values (computed as the difference of the 
double-precision and single-precision values divided by the double-precision 
value). In testl, as well as in all subsequent drivers, not all error flags are 
interrogated for possible malfunction. The user is urged, however, to do so 
as a matter of principle. 
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TABLE 3.1. Selected output from testl in the case 
of modified moments 

The routine 
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f mm (n , eps , modmom , om2 , f nu , ierr , f , f , rr ) 



used by the driver computes the first 2xn modified [ordinary] moments 
for up = om2, to a relative accuracy eps if modmom= . true . [.false.]. The 
results are stored in the array fnu. The arrays f , f 0, rr are internal working 
arrays of dimension n, and ierr is an error flag. On normal return, ierr = 
0; otherwise, ierr = 1, indicating lack of convergence (within a prescribed 
number of iterations) of the backward recurrence algorithm for computing 
the minimal solution {Cr(i^^)}. The latter is likely to occur iiuP is too close 
to 1. The routine fmm as well as its double-precision version dmm are listed 
immediately after the routine testl. 

In Table 3.1, we show selected results from the output of testl, when 
modmom= . true . (Complete results arc given in the package immediately 
after testl.) The values for A; = are expressible in terms of the complete 
elliptic integral, /?o = 2K(a;^), and were checked, where possible, against 
the 16S-values in Table 17.1 of [1]. In all cases, there was agreement to all 
16 digits. The largest relative error observed was 2.43 x 10~^^ for cj^ = .999 
and k = 2. When < .99^ the error was always less than 2.64 x 10 , 
which confirms the extreme stability of the modified Chebyshev algorithm 
in this example. It can be seen (as was to be expected) that for uP not too 
close to 1, the coefficients converge rapidly to ^. 

In contrast. Table 3.2 shows selected results (for complete results, see the 
package) in the case of ordinary moments (modmom= . false . ) and demon- 
strates the severe instability of the Chebyshev algorithm. We remark that 
the moments themselves are all accurate to essentially machine precision, as 
has been verified by additional computations. 



a;2 k 


err (3k 


u;2 


k 


err (3k 


.1 1 


1.187(-14) 


.9 


1 


3.270(-15) 


7 


2.603(-10) 




7 


4.819(-10) 


13 


9.663(-6) 




13 


1.841(-5) 


19 


4.251(-1) 




19 


6.272(-l) 


.5 1 


2.431 (-14) 


.999 


1 


6.311(-14) 


7 


5.571 (-10) 




7 


1.745 (-9) 


13 


9.307(-6) 




13 


8.589(~5) 


19 


5.798(-l) 




19 


4.808(0) 
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TABLE 3.2. Selected output from testl in the case of 
ordinary moments 

The next example deals with another weight function for which the mod- 
ified Chebyshev algorithm performs rather well. 

Example 3.2. dX^it) = t" In {l/t)dt on (0,1], a > -1. 
What is nice about this example is that both modified and ordinary 
moments of dX^ are known in closed form. The latter are obviously given 

by 

^kidX,) = , ! , , : fc= 0,1,2,..., (3.11) 

C7 + 1 + A; 

whereas the former, relative to shifted monic Legendre polynomials (ipoly=2 
in recur), are (cf. [17]) 

, ,,k-a <7^-'^{k-(J-l)\ 

(fc + ^ + 1)1 ' < a < k, a &n mteger, 

1 [ 1 \^(^L 1 \ ] A g + 1 - r 

otherwise. 

(3.12) 

The routines fmin and dmm appended to test2 in §1 of the package, simi- 
larly as the corresponding routines in Example 3.1, generate the first 2xn 
modified moments fQ, z^i, . . . , t'2n-i, if modmom= . true . , and the first 2xn 
ordinary moments, otherwise. The calling sequence of f mm is 

fmm(n,modmom,intexp,sigma,fnu) . 

The logical variable intexp is to be set .true., if a is an integer, and 
.false, otherwise. In either case, the input variable sigma is assumed of 
type real. 

The routine test 2 generates the first n recursion coefficients ak{dXu), 
Pk{dXcr) in single and double precision for a = — ^,0, |, where n = 100 
for the modified Chebyshev algorithm (modmom= . true . ), and n = 12 for the 
classical Chebyshev algorithm (modmom= . false . ). Selected double-precision 
results to 25 significant digits, when modified moments are used, are shown 
in Table 3.3. (The complete results are given in the package after test2.) 



VkidXa) = < 



15 



7 


k 


Oik 




.5 





.1111111111111111111111111 


4.000000000000000000000000 




12 


.4994971916094638566242202 


.06231277082877488477563886 




24 


.4998662912324218943801592 


.06245372557342242600457226 




48 


.4999652635485445800661969 


.06248855717748684742433618 




99 


.4999916184024356271670789 


.06249733823051821636937156 








.2500000000000000000000000 


1.000000000000000000000000 




12 


.4992831802157361310272625 


.06238356835953571123560330 




24 


.4998062839486146398501532 


.06247100084469111001639128 




48 


.4999494083797023879356424 


.06249281268110967462373889 




99 


.4999877992015903283047919 


.06249832670616925926204896 


.5 





.3600000000000000000000000 


.4444444444444444444444444 




12 


.4993755732917555644203267 


.06237082738280752611960887 




24 


.4998324497706394488722725 


.06246581011945496883543089 




48 


.4999567275223771727791521 


.06249115332711027176695932 




99 


.4999896931841789781887674 


.06249787251281682973825635 



TABLE 3.3. Selected output from test2 in the case of 
modified moments 

The largest relative errors observed, over all k = 0, 1, ... , 99, were respec- 
tively 6.211 x 10-l^ 2.237x10-^2^ 1.370x10-^2 for thea's, and 1.235x10-^°, 
4.446 X 10-^2^ 2.724 x IQ-^^ foj. ^j^g ^'g^ attained consistently at A; = 99. The 
accuracy achieved is slightly less than in Example 3.1, for reasons explained 
in [22, Example 5.3]. 

The complete results for a = — ^ are also available in [27, Appendix, 
Table 1]. (They differ occasionally by 1 unit in the last decimal place from 
those produced here, probably because of a slightly different computation of 
the modified moments.) The results for o" = can be checked up to /c = 15 
against the 30S-values given in [47, p. 92], and for 16 < < 19 against 
12S-values in [10, Table 3]. There is complete agreement to all 25 digits in 
the former case, and agreement to 12 digits in the latter, although there are 
occasional end figure discrepancies of 1 unit. These arc believed to be due 
to rounding errors committed in [10], since similar discrepancies occur also 
in the range k < 15. We do not know of any tables for a = ^, but a test 
will be given in §5, Example 5.1. 
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The use of ordinary moments (modinom= .false.) produces predictably 
worse results, the relative errors of which are shown in Table 3.4. 



k 


a err ak 


err Pk 


a err ak 


err Pk 


0" err 


err Pk 


2 


-.5 1.8(-13) 


7.7(-14) 


4.2(-13) 


7.6(-13) 


.5 1.6(-12) 


2.6(-13) 


5 


2.2(-9) 


1.2(-9) 


4.2(-9) 


1.2(-10) 


1.3(-8) 


6.6(-9) 


8 


l.l(-5) 


5.5(-6) 


4.3(-6) 


3.8(-6) 


6.0(-5) 


5.2(-6) 


11 


2.5(-l) 


1.7(-1) 


1.3(0) 


3.2(-l) 


2.2(0) 


4.7(-l) 



TABLE 3.4. Selected output from test2 in the case of 
ordinary moments 



4. STIELTJES, ORTHOGONAL REDUCTION, AND DISCRETIZATION PRO- 
CEDURES 

4.1. The Stieltjes procedure. It is well known that the coefficients ak{dX), 
Pk{d\) in the basic recurrence relation (1.3) can be expressed in terms of 
the orthogonal polynomials (1.2) and the inner product (1.1) as follows: 

ak{dX) = -, k>0- 

(7rfc,7rfe) ^^^^^ 

/3o(dA) = (7ro,7ro), Pk{dX) = . ^''"''''^ . , k>l. 

(7rfc_i,7rfc_ij 

Provided the inner product can be readily calculated, (4.1) suggests the 
following "bootstrap" procedure: Compute cto and /3o by the first relations 
in (4.1) for fc = 0. Then use the recurrence relation (1.3) for = to 
obtain vri. With ttq and tti known, apply (4.1) for = 1 to get ai, then 
again (1.3) to obtain 7r2, and so on. In this way, alternating between (4.1) 
and (1.3), we can bootstrap ourselves up to as many of the coefficients a^, 

as are desired. We attributed this procedure to Stieltjes, and called it 
Stieltjes's procedure in [19]. 

In the case of discrete orthogonal polynomials, i.e., for inner products of 
the form 

N 

{u,v) = ^Wku{xk)v{xk), Wk>0, (4.2) 
fe=i 
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Stieltjes's procedure is easily implemented; the resulting routine is called 
sti, and has the calling sequence 



St i (n , neap , X , w , alpha , bet a , ierr , pO , pi , p2) . 



On entry, 



neap 



x,w 



is the number of recursion coefficients desired; 
type integer 

is the number of terms, A'^, in the discrete in- 
ner product; type integer 

are arrays of dimension neap holding the ab- 
scissae X (k) = Xk and weights w (k) = Wk, k = 
1,2,..., neap, of the discrete inner product. 



On return. 



alpha, beta are arrays of dimension n containing the de- 
sired recursion coefficients alpha (k) = a^^i, 
beta {k) = /3k-i, A; = 1, 2, . . . , n 

ierr is an error flag having the value on normal 

return, and the value 1 if n is not in the proper 
range 1 < n < A'^; if during the computation 
of a recursion coefficient with index k there 
is impending underflow or overflow, ierr will 
have the value —A; in case of underflow, and 
the value +k in case of overflow. (No error 
flag is set in case of harmless underflow.) 

The arrays pO, pi, p2 arc working arrays of dimension neap. The double- 
precision routine has the name dsti . 

Occurrence of underflow [overflow] can be forestalled by multiplying all 
weights Wk by a sufficiently large [small] scaling factor prior to entering the 
routine. Upon return, the coefficient f3o will then have to be readjusted by 
dividing it by the same scaling factor. 
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4.2. Orthogonal reduction method. Another approach to producing the 
recursion coefficients a^, from the quantities Xk, Wk defining the inner 
product (4.2) is based on the observation (cf. [2], [29, §7]) that the symmetric 
tridiagonal matrix of order A'^ + 1, 







VP. 



'N-1 







aN-i 



(4.3) 



(the "extended Jacobi matrix" for the discrete measure dAjv impUed in 
(4.2)), is orthogonally similar to the matrix 



w 











" 





















XN . 



(4.4) 



Hence, the desired matrix J{d\N) can be obtained by applying Lanczos's 
algorithm to the matrix (4.4). This is implemented in the routine 



lancz (n , neap , x , w , alpha , beta , ierr , pO , pi ) , 



which uses a judiciously constructed sequence of Givens transformations 
to accomplish the orthogonal similarity transformation (cf. [45,3,41,2]; the 
routine lancz is adapted from the routine RKPW in [41, p. 328]). The 
input and output parameters of the routine lancz have the same meaning 
as in the routine sti, except that ierr can only have the value or 1, while 
pO, pi are again working arrays of dimension neap. The double-precision 
version of the routine is named dlancz. 

The routine lancz is generally superior to the routine sti: The pro- 
cedure used in sti may develop numerical instability from some point on, 
and therefore give inaccurate results for larger values of n. It furthermore is 
subject to underflow and overflow conditions. None of these shortcomings 
is shared by the routine lancz. On the other hand, there are cases where 
sti does better than lancz (cf. Example 4.5). 

We illustrate the phenomenon of instability (which is explained in [32]) 
in the case of the "discrete Chebyshev" polynomials. 
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Example 4.1. The inner product (4.2) with Xfc = — 1 + 2 jj-^, wj- = jj, 
k = l,2,...,N. 

This generates discrete analogues of the Legendre polynomials, which 
they indeed approach as N ^ oo. The recursion coefficients are explicitly 
known: 



ak = 0, k = 0,l,...,N -1; 

k = 1,2,...,N-1. 

To find out how well the routines sti and lancz generate them (in single 
precision), when A'^ = 40, 80, 160, 320, we wrote the driver testS, which 
computes the respective absolute errors for the a's and relative errors for 
the /3's. 

Selected results for Stieltjes's algorithm are shown in Table 4.1. The 
gradual deterioration, after some point (depending on A''), is clearly visible. 
Lanczos's method, in contrast, preserves essentially full accuracy; the largest 
error in the a's is 1.42(-13), 2.27(-13), 4.83(-13), 8.74(-13) for N = 40, 80, 
160, 320, respectively, and 3.38(-13), 6.63(-13), 2.17(-12), 5.76(-12) for the 
/?'s. 





n 


err a 


err f3 


A^ 


n 


err a 


err /3 


40 


< 35 


< 1.91(~13) 


< 7.78(-13) 


160 


< 76 


< 2.98(-13) 


< 7.61(-13) 




36 


3.01(-12) 


1.48(-11) 




85 


1.61 (-9) 


1.57(-8) 




37 


6.93(-ll) 


3.55(-10) 




94 


1.25(-4) 


1.17(-3) 




38 


2.57(-9) 


1.30 (-8) 




103 


2.64(-3) 


1.51(-1) 




39 


1.93(-7) 


9.58(-7) 




112 


2.35(-3) 


1.16(0) 


80 


< 53 


< 2.04(-13) 


< 6.92(-13) 


320 


< 106 


< 8.65(-13) 


< 7.39(-13) 




57 


2.04(-10) 


5.13("10) 




117 


3.96(-10) 


7.73("10) 




61 


3. 84 (-7) 


9.35(-7) 




128 


2.46(~6) 


4.67(-6) 




65 


1.94(-3) 


4.61(-3) 




139 


2.94(-2) 


6.27(-2) 




69 


1.87(-1) 


6.14(0) 




150 


1.15(-3) 


2.18(-2) 



TABLE 4.1. Errors in the recursion coefficients afc, Pk of (4.5) 
computed by Stieltjes's procedure 
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4.3. Multiple- component discretization procedure. We assume now a 
measure d\ of the form 

p 

dX{t) = w{t)dt + J2 ViKt - Xj)dt, P>0, (4.6) 
j=i 

consisting of a continuous part, w{t)dt, and (if p > 0) a discrete part written 
in terms of the Dirac (5-function. The support of the continuous part is 
assumed to be an interval or a finite union of disjoint intervals, some of 
which may extend to infinity. In the discrete part, the abscissae 
assumed pairwise distinct, and the weights positive, yj > 0. The inner 
product (1.1), therefore, has the form 




The basic idea of the discretization procedure is rather simple: One 
approximates the continuous part of the inner product, i.e., the integral in 
(4.7), by a sum, using a suitable quadrature scheme. If the latter involves 
terms, this replaces the inner product (4.7) by a discrete inner product 
(• ,-)n+p consisting oi N + p terms, the N "quadrature terms" and the 
p original terms. In effect, the measure dX in (4.6) is approximated by 
a discrete (N + p)-point measure dA^v+p- Wc then compute the desired 
recursion coefficients from the formula (4.1), in which the inner product 
(• , •) is replaced, throughout, by (• , •)n+p- Thus, in effect, we approximate 

ak{dX) ^ ak{dXN+p), Pk{dX) ^ l3k{dXN+p)- (4.8) 

The quantities on the right can be computed by the methods of §4.1 or §4.2, 
i.e., employing the routines sti or lancz. 

The difficult part of this approach is to find a discretization that results 
in rapid convergence, as N ^ oo, of the approximations on the right of (4.8) 
to the exact values on the left, even in cases where the weight function w 
in (4.6) exhibits singular behavior. (The speed of convergence, of course, is 
unaffected by the discrete part of the inner product (4.7).) To be successful 
in this endeavor often requires considerable inventiveness on the part of 
the user. Our routines, mcdis and dmcdis, that implement this idea in 
single rcsp. double precision, however, are designed to be flexible enough to 
promote the use of effective discretization procedures. 
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Indeed, if the support of the weight function w in (4.7) is contained in 
the (finite or infinite) interval (a, 6), it will often be useful to first decompose 
that interval into a finite number of subintervals, 



and to approximate the inner product separately on each subinterval [oj, bi], 
using an appropriate weighted quadrature rule. Thus, we write the integral 
in (4.7) as 



where Wi is an appropriate weight function on [aj,6j]. The intervals [aj,6j] 
are not necessarily disjoint. For example, the weight function w may be 
the sum w = wi + oi two weight functions on [a, b], which we may want 
to treat individually (cf. Example 4.2 below). In that case, one would 
take [ai,6i] = [027^2] = [0-,^] and wi on the first interval, W2 on the other. 
Alternatively, we may simply want to use a composite quadrature rule to 
approximate the integral, in which case (4.9) is a partition of [a, b] and 
Wi{t) = w{t) for each i. Still another example is a weight function w which 
is already supported on a union of disjoint intervals; in this case, (4.9) would 
be the same union, or possibly a refined union where some of the subintervals 
are further partitioned. 

In whichever way (4.9) and (4.10) are constructed, each integral on the 
right of (4.10) is now approximated by an appropriate quadrature rule, 



m 



supp w C [a,b] = [J [ai,bi], m > 1, 



(4.9) 



i=l 




(4.10) 




(4.11) 




(4.12) 



r=l 



This gives rise to the approximate inner product 




1=1 r=l 



(4.13) 



m 



1=1 
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In our routine mcdis, we have chosen, for simphcity, aU Ni to be the same 
integer Nq, 

Ni = No, z = l,2,...,m, (4.14) 

so that N = uiNq. Furthermore, if n is the number of ak and the num- 
ber of desired, we have used the following iterative procedure to de- 
termine the coefficients ak, Pk to a prescribed (relative) accuracy e: We 
let Nq increase through a sequence {iVQ'*'}s=o,i,2,... of integers, for each s 
use Stieltjes's (or Lanczos's) algorithm to compute = ak{dX^^^^[s]^^), 

f^k^ = (3k{dX [s] ), /c = 0, 1, . . . , n — 1, and stop the iteration for the first 
s > 1 for which all inequalities 

14'' - <el3^^\ = 0, 1, . . . , n - 1, (4.15) 

are satisfied. An error flag is provided if within a preset range Nq < iV™*^ 
the stopping criterion (4.15) cannot be satisfied. Note that the latter is 

based solely on the /3-coefRcients. This is because, unlike the a's, they are 
known to be always positive, so that it makes sense to insist on relative 
accuracy. (In our routine we actually replaced p]^^ on the right of (4.15) 
by its absolute value to insure proper termination in cases of sign-changing 
measures dX.) 

In view of the formulae (4.1), it is reasonable to expect, and indeed 
has been observed in practice, that satisfaction of (4.15) entails sufficient 
absolute accuracy for the a's if they are zero or small, and relative accuracy 
otherwise. 

Through a bit of experimentation, we have settled on the following se- 
quence of integers A^q*' : 

Ar['^l=2n, AW=iv[^-^l + A„ s = l,2,..., 

(4.16) 

Ai = l, A, = 2LV5J .n, s = 2,3,... . 

It should be noted that if the quadrature formula (4.11) is exact for 
each i, whenever u ■ v is a polynomial of degree < 2n — 1 (which is the 
maximum degree occurring in the inner products of (4.1), when k < n — 1), 
then our procedure converges after the very first iteration step! Therefore, 
if each quadrature rule Qi has (algebraic) degree of exactness > d{No), and 
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d{No)/No = 6 + 0{Nq^) as A^o ^ oo, then we let N^^^ = 1 + L(2n - l)/5\ 
in an attempt to get exact answers after one iteration. Normally, 5=1 (for 
interpolatory rules) or S = 2 (for Gauss-type rules). 

The calling sequence of the multiple-component discretization routine is 

mcdis (n , ncapm , mc , mp , xp , yp , quad , eps , iq, idelta , irout , 
f inl , f inr , endl , endr , xf er , wf er , alpha , bet a , neap , 
kount , ierr , ie , be , x , w , xm , wm , pO , pi , p2) 

dimension xp(*) ,yp(*) ,endl(mc) ,endr(mc) ,xf er(ncapm) , 
wfer(ncapm) ,alpha(n) ,beta(n) ,be(n) ,x(ncapm) , 
w(ncapm) ,xm(*) ,wm(*) ,pO(*) ,pl(*) ,p2(*) 

logical f inl, f inr 

On entry, 



n is the number of recursion coefficients desired; 

type integer 

ncapm is the integer Nq^^ above, i.e., the maximum 
integer A^o allowed (ncapm = 500 will usually 
be satisfactory) 

mc is the number of component intervals in the 
continuous part of the spectrum; type integer 

mp is the number of points in the discrete part of 
the spectrum; type integer; if the measure has 
no discrete part, set mp = 

xp , yp are arrays of dimension mp containing the ab- 
scissae and the jumps of the point spectrum 

quad is a subroutine determining the discretization 
of the inner product on each component in- 
terval, or a dummy routine if iq / 1 (see be- 
low); specifically, quad (n, x, w, i , ierr) pro- 
duces the abscissae x(r) = Xr,i and weights 
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w(r) = Wr,i, r = 1,2, ...,n, of the 77,-point 
discretization of the inner product on the in- 
terval [aj,6j] (cf. (4.13)); an error flag ierr is 
provided to signal the occurrence of an error 
condition in the quadrature process 



eps 



is the desired relative accuracy of the nonzero 
recursion coefficients; type real 



is an integer selecting a uscr-supplicd quadra- 
ture routine quad if iq = 1 or the ORTHPOL 
routine qgp (see below) otherwise 



idelta 



is a nonzero integer, typically 1 or 2, inducing 
fast convergence in the case of special quadra- 
ture routines; the default value is idelta = 1 



irout 



is an integer selecting the routine for generat- 
ing the recursion coefficients from the discrete 
inner product; specifically, irout = 1 selects 
the routine sti, and irout 7^ 1 the routine 
lancz. 



The logical variables f inl.f inr and the arrays endl.endr ,xf er ,wf er are 
input variables to the subroutine qgp and are used (and hence need to be 
properly dimensioned) only if iq 7^ 1. 

On return, 

alpha, beta are arrays of dimension n holding the desired 
recursion coefficients alpha (k) = ak-i, beta 
{k) = Pk-u k = 1,2, . . . , n 

neap is the integer A^o yielding convergence 

kount is the number of iterations required to achieve 

convergence 
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ierr is an error flag, equal to on normal return, 
equal to -1 if n is not in the proper range, 
equal to i if there is an error condition in the 
discretization on the ith interval, and equal 
to ncapm if the discretized Stieltjes procedure 
does not converge within the discretization 
resolution specified by ncapm 

ie is an error flag inherited from the routine sti 
or lancz (whichever is used). 

The arrays be ,x,w,xiii,wiii,pO ,pl ,p2 are used for working space, the last 
five having dimension mcxncapm+mp. 

A general-purpose quadrature routine, qgp, is provided for cases in which 
it may be difficult to develop special discretizations that take advantage of 
the structural properties of the weight function w at hand. The routine 
assumes the same setup (4.9)-(4.14) used in mcdis, with disjoint intervals 
[ai,bi], and provides for Qi in (4.12) the Fcjcr quadrature rule, suitably 
transformed to the interval [a^, bi], with the same number Ni = Nq of points 
for each i. Recall that the AT-point Fejer rule on the standard interval [-1,1] 
is the interpolatory quadrature formula 



where = cos((2r— l)7r/2A^) are the Chebyshev points. The weights are all 
positive and can be computed explicitly in terms of trigonometric functions 
(cf., e.g., [15]). The rule (4.17) is now apphed to the integral in (4.11) by 
transforming the interval [-1,1] to [ai,6i] via some monotone function (f)i (a 
linear function if [oj, 6^] is finite) and letting / = uvwf. 



N 



(4.17) 



r=l 




N 



r = l 



Thus, in effect, we take in (4.13) 



Xr,i = MXr)^ Wr,i = W^Wi{(t)i{x^))(t)'i{x^), Z = 1, 2, . . . , 



m. 



(4.18) 
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If the interval [0^,6^] is half-infinite, say of the form [0,cx)], we use (?!>i(t) = 
(1 + t)/{l — t), and similarly for intervals of the form [-00,6] and [a, 00]. If 
[oj, bi] = [—00, 00], we use 4>i{t) = t/{l — t^). 

The routine qgp has the following calling sequence: 

subroutine qgp (n , x , w , i , ierr , mc , f inl , f inr , endl , endr , xf er , wf er) 
dimension x(n) ,w(n) ,endl(mc) ,endr(mc) ,xf er(*) ,wf er(*) 
logical f inl, f inr 



On entry, 



is the number of terms in the Fejer quadrature 
rule 



i indexes the interval [ai,6i] for which the 

quadrature rule is desired; an interval that ex- 
tends to -00 has to be indexed by 1, one that 
extends to -l-oo by mc 

mc is the number of component intervals; type 
integer 

finl is a logical variable to be set .true, if the 
extreme left interval is finite, and .false, 
otherwise 



finr is a logical variable to be set .true, if the 
extreme right interval is finite, and .false, 
otherwise 



endl is an array of dimension mc containing the left 

endpoints of the component intervals; if the 
first of these extends to -00, endl(l) is not 
being used by the routine 

endr is an array of dimension mc containing the 
right endpoints of the component intervals; if 
the last of these extends to -|-oo, endr(mc) is 
not being used by the routine 
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xfer.wfer are working arrays holding the standard Fejer 
nodes and weights, respectively; they are di- 
mensioned in the routine mcdis. 

On return, 

x,w are arrays of dimension n holding the abscis- 

sae and weights (4.18) of the discretized inner 
product for the ith component interval 

ierr has the integer value 0. 

The routine calls on the subroutines fejer, symtr and tr, which are ap- 
pended to the routine qgp in §4 of the package. The first generates the Fejer 
quadrature rule, the others perform variable transformations. The user has 
to provide his own function routine wf (x, i) to calculate the weight function 
Wi{x) on the ith component interval. 

Example 4.2. Chebyshev weight plus a constant: w'^it) = {l — t^)~^^^ + c, 
c> 0, -1< t < 1. 

It would be difficult, here, to find a single quadrature rule for discretizing 
the inner product and obtain fast convergence. However, using in (4.9) 
m = 2, [ai,6i] = [02,62] = [-1,1], and wi{t) = (1 - t^)'^/^, W2{t) = c 
in (4.11), and taking for Qi the Gauss- Chebyshev, and for Q2 the Gauss- 
Legendre n-point rule (the latter multiplied by c), yields convergence to 
ak(w'^), I3k{w^), A: = 0, 1, . . . , n — 1, in one iteration (provided 5 is set equal 
to 2)! Actually, we need Nq = n + 1, m. order to test for convergence; cf. 
(4.15). The driver test4 implements this technique and calculates the first 
n = 80 beta-coefficients to a relative accuracy of 5000 xe* for c = 1, 10, 
100. (All aj- are zero.) Attached to the driver is the quadrature routine 
qchle used in this example. It, in turn, calls for the Gauss quadrature 
routine gauss, to be described later in §6. Anticipating convergence after 
one iteration, wc put ncapm = 81. 

The weight function of Example 4.2 provides a continuous link between 
the Chebyshev polynomials (c = 0) and the Legendre polynomials (c = 00); 
the recursion coefficients [3^ {w'^) indeed converge (except for A; = 0) to those 
of the Legendre polynomials, as c ^ 00. 

Selected results of test4 (where irout in mcdis can be arbitrary) are 
shown in Table 4.2. The output variable kount is 1 in each case, confirming 
convergence after one iteration. The coefficients l3o{w'^) are easily seen to be 
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TT + 2C. 



k 











5.141592654 


23.14159265 


203.1415927 


1 


.4351692451 


.3559592080 


.3359108398 


5 


.2510395775 


.2535184776 


.2528129500 


12 


.2500610870 


.2504824840 


.2505324193 


25 


.2500060034 


.2500682357 


.2501336338 


51 


.2500006590 


.2500082010 


.2500326887 


79 


.2500001724 


.2500021136 


.2500127264 



TABLE 4.2. Selected recursion coefficients Pkiw") for c = 1, 10, 100 



Example 4.3. Jacobi weight with one mass point at the left endpoint: 

u;("-/3)(^ -y) = [/if,"'^)]-i(l-x)°(l + x)^ + y 6{x + l) on (-1,1), ^uj,"'^) = 
2a+/3+ir(a + l)r(/3 + l)/r(a + /3 + 2), Q > -1, /3 > -1, y > 0. 

The recursion coefficients a^, Pk ^-re known explicitly (see [6, Eqs. (6.23), 
(3.5)]|^ and can be expressed, with some effort, in terms of the recursion 
coefficients aj^, for the Jacobi weight w^"'^\-) = (• ;0). The 

formulae are: 



Off — V 7 

1 + y 



_ J 2fc(a+fc) / _ n _L 2(/3+fc+l)(a+/3+fc+l) /J_ _ 

"A; — «fc -I- (Q,+^+2fc)(Q,+/3+2A;+l) '^''^ (a+/3+2A:+l)(a+/3+2A:+2) \ck 

(4.19) 



/?fc= ^4^, A; = 1,2,3,..., 



where 



1 + '-^^^^m^ydk 



co = l + y, Ck= — , A: = l,2,..., (4.20) 

1 + ydk 



^ In [6], the interval is taken to be [0,2], rather than [—1,1]. There is a typographical 
error in the first formula of (6.23), which should have the numerator 2/3 + 2 instead of 
2/3 + 1. 
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and 

di = 1, 

{l3 + k)ia + P + k) ^ 
= ^a + k-l){k-l) '^'-'^ ^ = ^' ^' 

Again, it is straightforward with mcdis to get exact results (modulo 
rounding) after one iteration, by using the Gauss- Jacobi quadrature rule 
(see gauss in §6) to discretize the continuous part of the measure. The 
driver test5 generates in this manner the first n = 40 recursion coefficients 
Ofc, Pk: k = 0,1, . . . ,n — 1, to a relative accuracy of SOOOxe*, for y = ^, 1, 2, 
4, 8. For each a = — .8(.2)1. and /? = -.8(.2)1., it computes the maximum 
relative errors (absolute error, if ~ 0) of the a^, Pk by comparing them 
with the exact coefficients. These have been computed in double precision 
by a straightforward implementation of the formulae (4.19)-(4.21). 

As expected, the output of test5 reveals convergence after one iteration, 
the variable kount having consistently the value 1. The maximum relative 
error in the Ofc is found to generally lie between 2 x 10^^ and 3 x 10^*^, the 
one m the Pk between 7.5 x 10"^^ and 8 x 10"^^; they are attained for k at 
or near 39. The discrepancy between the errors in the Ok and those in the 
Pk is due to the Ok being considerably smaller than the Pk, by 3-4 orders of 
magnitude. Replacing the routine sti in mcdis by lancz yields very much 
the same error picture. 

It is interesting to note that the addition of a second mass point at the 
other endpoint makes an analytic determination of the recursion coefficients 
intractable (cf. [6, p. 713]). Numerically, however, it makes no difference 
whether there are two or more mass points, and whether they are located in- 
side, or outside, or on the boundary of the support interval. It was observed, 
however, that if at least one mass point is located outside the interval [—1,1], 
the procedure sti used in mcdis becomes severely unstable]^ and must be 
replaced by lancz. 

Example 4.4. Logistic density function: w^x) = e~^/(l + e"^)"^ on 
(— oo, oo). 

^ This has also been observed in the similar Example 4.8 of [19], but was incorrectly 
attributed to a phenomenon of ill-conditioning. Indeed, the statement made at the end of 
Example 4.8 can now be retracted: stable methods do exist, namely the method embodied 
by the routine mcdis in combination with lancz. 
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(4.21) 



In this example we illustrate a slight variation of the discretization pro- 
cedure (4.9)-(4.13), which ends up with a discrete inner product of the same 
type as in (4.13) (and thus implementable by the routine mcdis) but derived 
in a somewhat different manner. The idea is to integrate functions with re- 
spect to the density w by splitting the integral into two parts, one from — oo 
to 0, the other from to oo, changing variables in the first part, and thus 
obtaining 

/oo roo p — t /'OO p — t 

JitMt)dt^l fi-t)^^^-^it + l (4.22) 

Since (1 + e~*)~^ quickly tends to 1 as t — > oo, a natural discretization of 
both integrals is provided by the Gauss-Laguerre quadrature rule applied to 
the product f{±t)- (1 -|- e~*)~^. This amounts to taking, in (4.13), m = 2 
and 

Xr,l = -Xr, Xr,2 = X^\ Wr,l = Wr,2 = ,^ , ^L.n > = 1, 2, . . . , A^, 

(1 -I- e~^r Y 

where x^, are the Gauss-Laguerre A^-point quadrature nodes and weights. 

The driver test 6 incorporates this discretization into the routines mcdis 
and dmcdis, runs them for n = 40 with error tolerances 5000xe* and 
1000 xe"', respectively, and prints the absolute errors in the a's (a^ = 0, 
in theory) and the relative errors in the /?'s. (We used the default value 
5 = 1.) Also printed are the number of iterations #it (= kount) in (4.15) 
and the corresponding final value (= neap). In single precision, we 
found #it = 1, A^Q =81, and in double precision, #it = 5, = 281. Both 
routines returned with the error flags equal to 0, indicating a normal course 



k 




err 


err Pk 





1 .000000000000000000000000 


4.572(-13) 


1.918(-13) 


1 


3.289868133696452872944830 


1.682(-13) 


5.641(-13) 


6 


89.44760352315950188817832 


2.187(-12) 


2.190(-12) 


15 


555.7827839879296775066697 


1.732(-13) 


2.915(-12) 


26 


1668.580222268668421827788 


3.772(-12) 


4.112(-12) 


39 


3753.534025194898387722354 


2.482(-ll) 


4.533(-12) 



TABLE 4.3. Selected output from teste 



31 



of events. A few selected double-precision values of the coefficients Pk along 
with absolute errors in the a's and relative errors in the /?'s are shown in 
Table 4.3. The results are essentially the same no matter whether sti or 
lancz is used in mcdis. The maximum errors observed are 2.482 x 10^^^ 
for the a's, and 4.939 x 10~^^ for the /3's, which are well within the single- 
precision tolerance e = 5000 x e^. 

On computers with limited exponent range, convergence difficulties may 
arise, both with sti and lancz, owing to underflow in many of the Laguerre 
quadrature weights. This seems to perturb the problem significantly enough 
to prevent the discretization procedure from converging. 

2 

Example 4.5. Half-range Hermite measure: 'w{x) = on (0, cxd). 

This is an example of a measure for which there do not seem to ex- 
ist natural discretizations other than those based on composite quadrature 
rules. Therefore, we applied our general-purpose routine qgp (and its double- 
precision companion dqgp), using, after some experimentation, the partition 
[0, oo] = [0, 3] U [3, 6] U [6, 9] U [9, oo]. The driver test? implements this, with 
n = 40 and an error tolerance 50 xe** in single precision, and 1000 xe'^ in 
double precision. 

The single-precision routine mcdis (using the default value 5=1) con- 
verged after one iteration, returning neap = 81, whereas the double-precision 



k 



Oik Pk 
err err fS^ 



.5641895835477562869480795 .8862269254527580136490837 
1.096(-13) 3.180(-13) 

1 .9884253928468002854870634 .1816901138162093284622325 
1.514(-13) 7.741(-14) 

6 2.080620336400833224817622 1.002347851011010842224538 

1.328(-13) 5.801(-14) 

15 3.214270636071128227448914 2.500927917133702669954321 

2.402(-14) 8.186(-14) 

26 4.203048578872001952660277 4.333867901229950443604430 

1.415(-13) 7.878(-14) 

39 5.131532886894296519319692 6.500356237707132938035155 

6.712(-13) 1.820(-14) 



TABLE 4.4. Selected output from test7 
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routine dmcdis took four iterations to converge, and returned ncapd = 201. 
Selected results (where err and err (3^ both denote relative errors) are 
shown in Table 4.4. The maximum error err occurred at A; = 10, and 
had the value 1.038 x 10~^^, whereas max^ err I3k = 3.180 X 10-^3 is at- 
tained at A; = 0. The latter is within the error tolerance e, the former only 
slightly larger. Comparison of the double-precision results with Table I on 
the microfiche supplement to [12] revealed agreement to all 20 decimal dig- 
its given there, for all k in the range < /c < 19. Interestingly, the routine 
sti in mcdis did consistently better than lancz on the /3's, by a factor as 
large as 235 (for k = 33), and is comparable with lancz (sometimes better, 
sometimes worse) on the a's. 

Without composition, i.e., using mc=l in mcdis, it takes 8 iterations 
(Nq = 521) in single precision, and 10 iterations {Nq = 761) in double 
precision, to satisfy the much weaker error tolerances £ = 5 10~^ and e*^ = 
^ 10~^^, respectively. All single-precision results, however, turn out to be 
accurate to about 12 decimal places. (This is because of the relatively large 
final increment Ag = 2n = 80 in Nq (cf. (4.16)) that forces convergence.) 

4.4. Discretized modified Chebyshev algorithm. The whole apparatus of 
discretization (cf. (4.9)-(4.14)) can also be employed in connection with the 
modified Chebyshev algorithm (cf. §3), if one discretizes modified moments 
rather than inner products. Thus, one approximates (cf. (4.14), (4.16)) 

Uk{d\) « Md)^^j,ls]^p) (4.23) 

and iterates the modified Chebyshev algorithm with s = 0,1,2,... until 

the convergence criterion (4.15) is satisfied. (It would be unwise to test 
convergence on the modified moments, for reasons explained in [19, §2.5].) 
This is implemented in the routine mccheb, whose calling sequence is 

mccheb (n , ncapm , mc , mp , xp , yp , quad , eps , iq , idelt a , f inl , 
f inr , endl , endr , xf er , wf er , a , b , f nu , alpha , bet a , neap , 
kount , ierr ,be,x,w,xiii,wiii,s,s0,sl,s2) 

Its input and output parameters have the same meaning as in the routine 
mcdis. In addition, the arrays a,b of dimension 2xn-l are to be supplied 
with the recursion coefficients a.{k) = o,fc-i, b(A;) = bk-i, k = 1,2, . . . , 2xn- 
1, defining the modified moments. The arrays be,x,w,xm,wm,s,s0,sl,s2 
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are used for working space. The double-precision version of the routine has 
the name dmcheb. 

The discretized modified Chebyshev algorithm must be expected to be- 
have similarly as its close relative, the modified Chebyshev algorithm. In 
particular, if the latter suffers from ill-conditioning, so does the former. 

Example 4.6. Example 3.1, revisited. 

We recompute the n = 40 first recursion coefficients a^, Pk of Example 
3.1 to an accuracy of lOOxe* in single precision, using the routine mccheb 
instead of the routine cheb. For the discretization of the modified moments 
we employed the Gauss-Chebyshev quadrature rule: 

I m{i-uHY'/\i-tY'/'dt^ ^ Y.fi^r){i-^'xi)-y\ (4.24) 

1 r=l 

where = cos((2r — 1)tt/2N) are the Chebyshev points. This is accom- 
plished by the driver testS. The results of this test (shown in the package) 
agree to all 10 decimal places with those of testi. The routine mccheb 
converged in one iteration, with neap = 81, for = .1, .3, .5, .7, .9, in 4 
iterations, with neap = 201, for cj^ = .99, and in 8 iterations, with neap = 
521, for a;^ = .999. A double-precision version of tests was also run with 
£ = ^ X 10~^° (not shown in the package) and produced correct results to 
20 decimals in one iteration (neap = 81) for = .1, .3, .5, .7, in 3 iterations 
(neap = 161) for u;^ = .9, in 6 iterations (neap = 361) for uP' = .99, and in 

II iterations (neap = 921) for oj"^ = .999. 



5. MODIFICATION ALGORITHMS 

Given a positive measure d\{t) supported on the real line, and two poly- 
nomials u{t) = lb Wp^^{t — Up), v{t) = n^^^(i — Va) whose ratio is finite on 
the support of d\, we may ask for the recursion coefficients = a!jt(dA), 
(^k = {d^) of the modified measure 

d\{t) = ^ dX{t), t e supp(dA), (5.1) 

v[t) 

assuming known the recursion coefficients ak = ak{dX), Pk = PkidX) of the 
given measure. Methods that accomplish the passage from the a's and /?'s 
to the a's and /3's will be called modification algorithms. The simplest case 
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s = (i.e., v{t) = 1) and u positive on supp(ciA) has already been considered 
by ChristofFel [7], who represented the polynomial n(-)7rfc(-) = n(-)7rfc (• ; dX) 
in determinantal form in terms of the polynomials vrj(-) = Trj{- ;dX), j = 
k,k + 1, . . . ,k + r. This is now known as Christoffel's theorem,. Christoffel, 
however, did not address the problem of how to generate the new coefficients 
ak, $k in terms of the old ones. For the more general modification (5.1), 
Christoffel's theorem has been generalized by Uvarov [48,49]. The coefhcient 
problem stated above, in this general case, has been treated in [20], and 
previously by Galant [13] in the special case v{t) = 1. 

The passage from dX to dX can be carried out in a sequence of elementary 
steps involving real linear factors t — x,OT real quadratic factors {t — x)'^ + y'^, 
either in u{t) or in v{t). The corresponding elementary steps in the passage 
from the a's and /?'s to the d's and /3's can all be performed by means of 
certain nonlinear recurrences. Some of these, however, when divisions of the 
measure dX are involved, are liable to instabilities. An alternative method 
can then be used, which appeals to the modified Chebyshev algorithm sup- 
plied with appropriate modified moments. These latter are of independent 
interest and find application, e.g., in evaluating the kernel in the contour 
integral representation of the Gauss quadrature remainder term. 

5.1 Nonlinear recurrence algorithms. The routine that carries out the 
elementary modification steps is called chri and has the calling sequence 

chri (n , iopt , a , b , x , y , hr , hi , alpha , beta , ierr ) . 

On entry, 

n is the number of recursion coefficients desired; 

type integer 
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iopt is an integer identifying the type of modifica- 
tion as follows: 
1: dX{t) = (t - x)dXit) 
2: dX{t) = i{t-xf + y^)d\{t),y>0 
3: dX = + y'^)dX{t) with dX{t) and 
supp(dA) assumed symmetric with 
respect to the origin and y > 
4: dX{t) = dX{t)/{t - x) 
5: dX{t) = dX{t)/{{t-xf + y'^),y>Q 
6: dX{t) = dX{t)/{f + y"^) with dX{t) and 
supp(dA) assumed symmetric with 
respect to the origin and y > 
7: dX{t) = {t- xfdX{t) 
a,b are arrays of dimension n+1 holding the re- 
cursion coefficients a.{k) = a;fc_i(dA), b(A;) = 
h-i{dX), A; = 1,2,..., n+1 
x,y are real parameters defining the linear and 
quadratic factors (or divisors) of dX 



hr.hi are the real and imaginary part, respectively, 
of /j^(iA(t)/(z — t), where z = x + iy; the 
parameter hr is used only if iopt = 4 or 5, 
the parameter hi only if iopt = 5 or 6. 



On return, 



alpha, beta are arrays of dimension n containing the de- 
sired recursion coefficients alpha(/c) = ak-i 
(dX), heta{k) = Pk-i (dX), = 1, 2, . . . , n 

ierr is an error flag, equal to on normal return, 

equal to 1 if n < 1 (the routine assumes that 
n is larger than or equal to 2), and equal to 2 
if the integer iopt is inadmissible. 

It should be noted that in the cases iopt = 1 and iopt = 4, the modified 
measure dX is positive [negative] definite if x is to the left [right] of the 
support of dX, but indefinite otherwise. Nevertheless, it is permissible to 
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have X inside the support of dX (or inside its convex hull), provided the 
resulting measure d\ is still quasi-definite (cf. [20]). 

For iopt = 1,2, ...,6, the methods used in chri are straightforward 
implementations of the nonlinear recurrence algorithms respectively in Eq. 
(3.7), (4.7), (4.8), (5.1), (5.8) and (5.9) of [20]. The only minor modification 
required concerns Po = Po{dX). In [20], this constant was taken to be 0, 
whereas here it is defined to be $o = /j^dA(t). Thus, for example, if iopt 
= 2, 

Po= f {{t - x f + y^)dX{t) = f ({t-ao + ao-x f + y^)dX{t) 



{{t - ao) + (ao - xY + i/)dX{t), 

R. 

since /j^(t — ao)dA(t) = /j^ 7ri(i)dA(t) = 0. Furthermore (cf. (4.1)), 

/ {t-aofdX{t)=Mi, 
so that the formula to be used for (3q is 

/3o = Wi + («o - xf + y^) (iopt = 2) . 

Similar calculations need to be made in the other cases. 

The case iopt = 7 incorporates a QR step with shift x, following Kautsky 
and Golub [42], and uses an adaptation of the algorithm in [51, Eq. (67.11), 
p. 567] to carry out the QR step. The most significant modification made 
is the replacement of the test c 7^ by |c| > e, where £ = 5 x is a 
quantity close to, but slightly larger than, the machine precision. (Without 
this modification, the algorithm could fail.) 

The methods used in chri are believed to be quite stable when the 
measure dX is modified multiplicatively (iopt = 1, 2, 3 and 7). When 
divisions are involved (iopt = 4, 5 and 6), however, the algorithms rapidly 
become unstable as the point z = x + iy G C moves away from the support 
interval of dX. (The reason for this instability is not well understood at 
present; see, however, Galant [14].) For such cases there is an alternative 
routine, gchri (see §5.2), that can be used. 
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Example 5.1. Checking the results (for o" = 2) of test2. 

We apply chri (and the corresponding double precision routine dchri) 
with iopt = 1, X = 0, to d\a{t) = In {l/t) on (0,1) with a = -\, to 
recompute the results of test2 for cr = ^. This can be done by a minor 



.5 



k 


12 
24 
48 
98 



err ak 

7.895(-14) 
3.280(-12) 
7.648(-12) 
2.076(-ll) 
6.042(-ll) 



err f3k 

4.796(-14) 
6.195(-12) 
1.478(~11) 
4.088(-ll) 
1.201(-10) 



err 

2.805(-28) 
8.958(-26) 
2.065(-25) 
5.683(-25) 
1.504(-24) 



err /3j 



7.952(-28) 
1.731(-25) 
3.985(-25) 
1.121(-24) 
2.987(-24) 



TABLE 5.1. Comparison between modified Chebyshev algorithm 
and modification algorithm in Example 5.1 (cf. Example 3.2) 



modification, named test9, of test2. Selected results from it, showing 
the relative discrepancies between the single-precision values a^, resp. 
double-precision values a^, computed by the modified Chebyshev algo- 
rithm and the modification algorithm, are shown in Table 5.1 (cf. Table 
3.3). The maximum errors occur consistently for the last value of k (= 98). 

Example 5.2. Induced orthogonal polynomials. 

Given an orthogonal polynomial iTmi' ]dX) of fixed degree m > 1, the 
sequence of orthogonal polynomials %,,„(•) = 7rfc(- ;7r^dA), A; = 0, 1, 2, . . . , 
has been termed induced orthogonal polynomials in [33]. Since their measure 
dXm modifies the measure dX by a product of quadratic factors, 

m 

dXmit) = l[it-x^f-dXit), (5.2) 

where are the zeros of tt^, we can apply the routine chri (with iopt=7) 
m times to compute the n recursion coefficients ak,m = ctkidXm), $k,m = 
/3k{dXrn), k = 0,1, . . . ,n — 1, from the n + m coefficients = ak{dX), (3k = 
(3k{dX), k = 0,1, . . . , n—l+m. The subroutines indp and dindp in the driver 
testlO implement this procedure in single resp. double precision. The 
driver itself uses them to compute the first n = 20 recursion coefficients of the 
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induced Legendre polynomials with m = 0, 1, . . . , 11. It also computes the 
maximum absolute errors in the a's {ak,m = for all m) and the maximum 
relative errors in the ^'s by comparing single-precision with double-precision 
results. 

An excerpt of the output of test 10 is shown in Table 5.2. It already 
suggests a high degree of stability of the procedure employed by indp. This 
is reinforced by an additional test (not shown in the package) generating n 
= 320 recursion coefficients (3k,nn < /c < 319, for m = 40, 80, 160, 

320 and dX being the Legendre, the Ist-kind Chebyshev, the Laguerre, and 
the Hermite measure. Table 5.3 shows the maximum absolute error in the 
Q!fe,m) < ^ < 319 (relative error in the Laguerre case), and the maximum 
relative error in the $k,mj < k < 319. 



k 


m = 0, Pk,m 


m = 2, Pk,m 


m = 6, Pk,m 


m = 11, Pk,m 





2.0000000000 


.1777777778 


.0007380787 


.0000007329 


1 


.3333333333 


.5238095238 


.5030303030 


.5009523810 


6 


.2517482517 


.1650550769 


.2947959861 


.2509913424 


12 


.2504347826 


.2467060415 


.2521022519 


.1111727541 


19 


.2501732502 


.2214990335 


.2274818789 


.2509466619 



err a 0.000(0) 
err P 1.737(~14) 



1.350(~13) 
2.032(-13) 



9.450(-13) 
2.055(-12) 



1.357(-12) 
3.748(-12) 



TABLE 5.2. Induced Legendre polynomials 



Legendre Chebyshev Laguerre Hermite 

m err a err $ err a err $ err a err $ err a err $ 

40 3.4(-ll) 1.5(-10) 1.9(-9) 7.9(-10) 3.0(-10) 6.0(-10) 1.8(-9) 2.7(-10) 

80 1.4(-10) 5.4(-10) 2.1(-9) 2.2(-9) 5.8(-10) 9.2(-10) 7.9(-9) 9.2(-10) 

160 1.5(-9) 5.1(-9) 9.5(-9) l.l(-8) 7.8(-10) 1.4(-9) l.l(-8) 6.8(-10) 

320 3.3(-9) 2.1(-8) 9.6(-9) 2.1(-8) 3.9(-9) 7.2(-9) 2.5(-8) l.l(-9) 



TABLE 5.3. Accuracy of the recursion coefficients for 
some classical induced polynomials 
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5.2. Methods based on the modified Chebyshev algorithm. As was noted 
earlier, the procedure chri becomes unstable for modified measures involv- 
ing division of dX{t) by t — x or (t — x)^ +y'^asz = x + iyGC moves away 
from the "support interval" of dX, i.e., from the smallest interval contain- 
ing the support of dX. We now develop a procedure that works better the 
further away z is from that interval. 

The idea is to use modified moments of dX relative to the polynomi- 
als 7rfc(- ;dX) to generate the desired recursion coefficients a^, (3^ via the 
modified Chebyshev algorithm (cf. §3). The modified moments in question 
are 

Vk = Uk{x-dX)= [ ^^^^2^dX{t), A; = 0,1,2,..., (5.3) 
for linear divisors, and 

Uk = i^kix,y;dX)= [ dX{t), A; = 0, 1,2, . . . , (5.4) 

for quadratic divisors. Both can be expressed in terms of the integrals 

pj, = Pk{z;dX)= [ ^feM^) dX{t), C\supp(ciA), fc = 0,1,2,..., 

(5.5) 

the first by means of 

Vk{x]dX) = -pk{z\dX), z = x, (5.6) 
the others by means of 

, , Im pdz] dX) . ,^ „, 

ffc(x,y;dA) = , z = x + iy. (5.7) 

Im z 

The point to observe is that {^Ph{z\ dX)\ is a minimal solution of the basic 
recurrence relation (1.3) for the orthogonal polynomials {vrfe(- \dX)) (cf. 
[18]). The quantities pk{z;dX), k = 0, 1, . . . ,n, therefore can be computed 
accurately by a backward recurrence algorithm ([18, §5]) which, for u > n, 
produces approximations p^l^\z;dX) converging to pk{z;dX) when i/ — oo, 
for any fixed k, 

p^^\z;dX) pkiz;dX), v ^ oo. (5.8) 
The procedure is implemented in the routine 



knum (n , nuO , numax , z , eps , a , b , rho , nu , ierr , rold) , 
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which computes Pkiz] dX) for A; = 0, 1, . . . , n to a relative precision eps. The 
results are stored as rho(/c) = pk_i{z; dX), k = 1,2, . . . , n+1, in the complex 
array rho of dimension n + 1. The user has to provide a starting index nuO 
= > n for the backward recursion, which the routine then increments 
by units of 5 until convergence to within eps is achieved. If the requested 
accuracy eps cannot be realized for some v < numax, the routine exits with 
ierr = numax. Likewise, if uq > numax, the routine exits immediately, with 
the error flag ierr set equal to nuO. Otherwise, the value of v for which 
convergence is obtained is returned as output variable nu. The arrays a, b 
of dimension numax are to hold the recursion coefficients a(A;) = a^-iidX), 
h{k) = Pk-i{dX), k = 1,2,..., numax, for the given measure dX. The 
complex array rold of dimension n + 1 is used for working space. In the 
interest of rapid convergence, the routine should be provided with a realistic 
estimate of uq. For classical measures, such estimates are known (cf. [18, 
§5]) and are implemented here by the function routines 



nuOjac(n,z,eps) , nu01ag(n,z,al,eps) , nuOher(n,z,eps). 



The first is for Jacobi measures, the second for generalized Laguerre mea- 
sures with parameter al = a > — 1, and the last for the Hermite measure. 
Note that uq for Jacobi measures does not depend on the weight parameters 
a, (3, in contrast to vq for the generalized Laguerre measure. 

The name knum comes from the fact that Pn{z', dX) in (5.5) is the numer- 
ator in the kernel 

/ ^^^ Pn{z;dX) 
Kn{z;dX) = , . (5.9) 

TTn[z; dX) 

of the remainder term of the n-point Gaussian quadrature rule for analytic 
functions (cf., e.g., [36]). For the sequence of kernels Kq, Ki,...,Kn we 
have the following routine: 



subroutine kern (n , nuO , numax , z , eps , a, b , ker ,nu , ierr , rold) 

complex z, ker , rold, pO,p,pml 

dimension a(numax) ,b(numax) ,ker(*) ,rold(*) 

call knum (n , nuO , numax , z , eps , a , b , ker , nu , ierr , rold) 

if (ierr .ne . 0) return 

pO=(0. ,0.) 

p=(l. ,0.) 

do 10 k=l,n 
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pnil=pO 
pO=p 

p= (z-a (k) ) *pO-b (k) *pml 
ker(k+l)=ker(k+l)/p 
10 continue 
return 
end 

The meaning of the input and output parameters is the same as in knum. 
The double precision version of the routine is named dkern. 

All the ingredients are now in place to describe the workings of gchri, 
the alternative routine to chri when the latter is unstable. First, the rou- 
tine knum is used to produce the first 2n modified moments Uk{x;dX) resp. 
VkixjU; dX), k = 0, 1, . . . , 2n — 1. These are then supplied to the routine 
cheb along with the recursion coefficients ak{dX), Pk{dX) (needed anyhow for 
the computation of the u^), which produces the desired coefficients ak{dX), 
Pk{dX), k = 0,1, . . . ,n — 1. The routine has the calling sequence 

gchr i (n , i opt , nuO , numax ,eps,a,b,x,y, alpha , bet a , 
nu , ierr , ierr c , f nu , rho ,rold,s,s0,sl,s2) . 

On entry, 

n is the number of recursion coefficients desired; 

type integer 

iopt is an integer identifying the type of modifica- 
tion as follows: 

1: dX{t) = dX{t)/{t — x), where x is assumed 
outside the smallest interval containing 
supp((iA) 

2: dX{t) = dX{t)/{{t-xf + y'^),y>Q 

nuO is an integer vq > 2n estimating the starting 
index for the backward recursion to compute 
the modified moments; if no other choices are 
available, take nuO = 3xn 
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numax 



is an integer used to terminate backward re- 
cursion in case of nonconvergence; a conser- 
vative choice is nvmiax=500 



eps is a relative error tolerance; type real 

a,b are arrays of dimension numax to be sup- 

plied with the recursion coefficients a{k) = 
ak-i{d\), b(fc) = h-i{dX), k = 1,2,..., 
numax 

x,y are real parameters defining the linear and 

quadratic divisors of dX. 

On return, 

alpha, beta are arrays of dimension n containing the de- 
sired recursion coefficients alpha(A;) = ak-i 
beta(A;) = k = 1,2, . . . , n 

nu is the index for which the error tolerance 

eps is satisfied for the first time; if it is never 
satisfied, nu will have the value numax 

ierr is an error flag, where 

ierr = on normal return 

ierr = 1 if iopt is inadmissible 

ierr = nuO if nuO > numax 

ierr = numax if the backward recurrence 
algorithm does not converge. 

ierr = -1 if n is not in range 
ierrc is an error flag inherited from the routine 

cheb. 

The real arrays fnu,s,s0,sl,s2 are working space, all of dimension 2xn 
except s, which has dimension n. The complex arrays rho, rold are also 
working space, both of dimension 2n. The routine calls on the subroutines 
knum and cheb. The double-precision version of gchri has the name dgchri. 
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Since the routine gchri is based on the modified Chebyshev algorithm, 
it shares with the latter its proneness to ill-conditioning, particularly in 
cases of measures supported on an infinite interval. On finitely supported 
measures, however, it can be quite effective, as will be seen in the next 
example. 

Example 5.3. The performance of chri and gchri. 

To illustrate the severe limitations of the routine chri in situations where 
divisions of the measure dX are involved, and at the same time to document 
the effectiveness of gchri, we ran both routines with n = 40 for Jacobi 
measures dX^°''^^ with parameters a, /3 = — .8(.4).8, P > a. This is done in 
testll. 

The routine testll first tests division by t — x, where x = 1.001, 
-1.01 , -1.04, -1.07, -1.1. Both routines chri and gchri are run in sin- 
gle and double precision, the latter with e = 10 x and e = 100 x e'^, 
respectively. The double-precision results are used to determine the abso- 
lute errors in the a's and the relative errors in the /3's for each routine. The 
required coefficients a/c, /3a;, < A; < fmax — 1 (j^max = 500 for single preci- 
sion, and 800 for double precision) are supplied by recur and drecur with 
ipoly = 6. The routine nuOjac is used to provide the starting recurrence 
index uq resp. fg. In Tables 5.4 and 5.5, relating respectively to linear 
and quadratic divisors, we give only the results for the Legendre measure 
{a = P = 0). The first line in each 3- line block of Table 5.4 shows x, uq, Uq 
and the maximum (over k, < k < 39) errors in the and Pk for gchri, 
followed by the analogous information (excepting the i^o's) for chri. The 
recurrence index v yielding convergence was found (not shown in testll) 
to be 1/ = i/Q -|- 5 and = + 5, without exception. 
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gchri chri 
X uq Vq err a err /3 err a err /3 



-1.001 


418 757 
reconstruction 
errors 


8.000( 

8.527( 
1.421( 


-14) 

-14) 
-14) 


1.559( 

1.705( 
5.329( 


-13) 

-13) 
-14) 


1.013( 

1.010( 
2.019( 


-13) 

-27) 
-28) 


1.647( 

2.423( 
1.211( 


-13) 

-27) 
-27) 


-1.01 


187 


294 


4.016( 
3.553( 
7.105( 


-14) 
-14) 
-15) 


6.907( 
9.946( 
4.262( 


-14) 
14) 
-14) 


1.396( 

6.058( 
1.515( 


-10) 

-28) 
-28) 


2.424( 
1.211( 
9.080( 


-10) 

-27) 
-28) 


-1.04 


133 


187 


3.590( 
2.842( 
7.105( 


-14) 
-14) 
-15) 


4.759( 
7.103( 
4.263( 


-14) 
-14) 
-14) 


5.944( 
5.554( 
1.010( 


-6) 

-28) 
-28) 


8.970( 

1.312( 
9.080( 


-6) 

-27) 
-28) 


-1.07 


120 


161 


2.194( 
2.842( 
7.105( 


-14) 
-14) 
-15) 


4.850( 
7.104( 
4.263( 


-14) 
-14) 
-14) 


5.334( 
6.058( 
1.010( 


-3) 

-28) 
-28) 


7.460( 
1.211( 
7.062( 


-3) 

-27) 

-28) 


-1.1 


114 


148 


2.238( 
2.132( 
1.549( 


-14) 
-14) 
-12) 


4.359( 
5.683( 
1.833( 


-14) 
-14) 
-12) 


4.163( 
3.534( 
1.010( 


0) 

-28) 
-28) 


4.959( 
1.009( 
6.057( 


+1) 

-27) 

-28) 



TABLE 5.4. Performance of gchri and chri for elementary 
divisors t — x oi the Legendre measure d\{t) 



It can be seen from the leading lines in Tabic 5.4 that chri rapidly loses 
accuracy as x moves away from the interval [—1,1], all single-precision ac- 
curacy being gone by the time x reaches -1.1. Similar, if not more rapid, 
erosion of accuracy is observed for the other parameter values of a, (3. The 
next two lines in each 3-line block show "reconstruction errors", i.e., the 
maximum errors in the ct's and /3's if the a's and /3's produced by gchri, 
chri and dgchri, dchri are fed back to the routines chri and dchri with 
iopt = 1 to recover the original recursion coefficients in single and double 
precision. The first of these two lines shows the errors in reconstructing these 
coefficients from the output of gchri (resp. dgchri), the second from the 
output of chri (resp. dchri). Rather remarkably, the coefficients are recov- 
ered to essentially full accuracy, even when the input coefficients (produced 
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by chri and dchri) are very inaccurate! This is certainly a phenomenon 
that deserves further study. It can also be seen from Table 5.4 (and the 
more complete results in §1 of the package) that gchri consistently pro- 
duces accurate results, some shght deterioration occurring only very close 
to a; = — 1, where the routine has to work harder. 

The second half of testll tests division by (t — x)^ + ?/^, where z = x + iy 
is taken along the upper half of the ellipse 

£p = {zeC: z = ^ (^pe'"^ + ^e"'^^ , 0<^< 27r}, p>l, (5.10) 

which has foci ±1 and sum of the semiaxes equal to p. (These ellipses are 
contours of constant uq for Jacobi measures.) We generated information 
analogous to the one in Table 5.4, for p = 1.05, 1.1625, 1.275, 1.3875, 1.5, 
except that all quantities are averaged over 19 equally spaced points on £p 
corresponding to = j7r/20, j = 1, 2, . . . , 19. Selected results (bars indicate 
averaging), again for the Legendre case, are shown in Table 5.5. They reveal 
a behavior very similar to the one in Table 5.4 for linear divisors. 











gch: 


ri 




chri 


p 






err 


a 


eff 




err a 


err (5 


1.05 


390 700 

reconstruction 
errors 


7.879(- 
7.814(^ 
2.024(- 


-13) 

"13) 
-14) 


1.440(- 
1.433(- 
8.442(- 


-12) 
-12) 
-14) 


7.685(-14) 
1.768(-26) 
3.016(-28) 


1.556(-13) 
3.042(-26) 
1.742 (-27) 


1.275 


142 


204 


6.252(- 
6.554(- 
2.295(- 


-14) 
-14) 
-14) 


1.287(- 
1.279(- 
8.970(- 


-13) 
-13) 
-14) 


4.562(-7) 

1.541(-27) 

3.579(-28) 


6.162(-7) 
3.061( 27) 
1.646(-27) 


1.5 


117 


154 


3.991( 

4.207(- 
3.805(- 


-14) 
-14) 
-14) 


7.966(- 
9.064(- 
8.971(- 


-14) 
-14) 
-14) 


4.906(-l) 

6.932(-28) 

4.351(-28) 


2.339(0) 

1.676(-27) 

1.744(-27) 



TABLE 5.5. Performance of gchri and chri for elementary 
divisors {t — x)"^ + of the Legendre measure d\{t) with z = x + iy on £p. 
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6. GAUSS-TYPE QUADRATURE RULES 



One of the important uses of orthogonal polynomials is in the approxi- 
mation of integrals involving a positive measure d\ by quadrature rules of 
maximum, or nearly maximum, algebraic degree of exactness. In this con- 
text, it is indispensable to know the recursion coefficients for the respective 
orthogonal polynomials {7rfc(- ; dX)}, since they allow us to generate the de- 
sired quadrature rules accurately and effectively via eigenvalue techniques. 
The software developed in the previous sections thus finds here a vast area 
of application. 



6.1. Gaussian quadrature. Given the (positive) measure d\ (having an 
infinite number of support points), there exists, for each n G N, a quadrature 
rule 



R 



f{t)d\{t) = Y,Wkf{xu) + R^{f) 



(6.1) 



fe=i 



having algebraic degree of exactness 2n — 1, i.e., zero error, Rn{f) = 0, 
whenever / is a polynomial of degree < 2n — 1. The nodes Xk indeed are the 
zeros of the nth-degree orthogonal polynomial 7r„(- ;c?A), and the weights 
Wk, which are all positive, are also expressible in terms of the same orthog- 
onal polynomials. Alternatively, and more importantly for computational 
purposes, the nodes Xk are the eigenvalues of the nth-order Jacobi matrix 



JnidX) 



"0 







^2 







\//3n-l "n-l 



(6.2) 



where ol^ = Ofc(f^A), = PkidX) are the recurrence coefficients for the 
(monic) orthogonal polynomials {7rfc(- ]dX)}, and the weights Wk are ex- 
pressible in terms of the associated eigenvectors. Specifically, if 



Jn{d\)vk = XkVk, v^Vk = 1, = 1, 2, . . . , n, 



(6.3) 



i.e., Vk is the normalized eigenvector of Jn{dX) corresponding to the eigen- 
value Xk, then 

Wk = l3ovl^i, k = l,2,...,n, (6.4) 



47 



where /3o = Poid^) is defined in (1.4) and Vk,i is the first component of 
Vk (cf. [40]). There are well-known and efficient algorithms, such as the 
QR algorithm, to compute eigenvalues and (part of the) eigenvectors of 
symmetric tridiagonal matrices. These are used in the routine gauss^, whose 
calling sequence is 



gauss (n , alpha , beta , eps , zero , weight , ierr , e) . 



On entry, 

n is the number of terms in the Gauss formula; 

type integer 

alpha, bet a are arrays of dimension n assumed to hold 
the recursion coefficients alpha(A;) = ak-i, 
beta(/c) = Pk-i, k = 1,2, . . . , n 

eps is a relative error tolerance, for example, the 

machine precision. 

On return, 

zero , weight are arrays of dimension n containing the nodes 
(in increasing order) and the corresponding 
weights of the Gauss formula, zero (A;) = Xk, 
weight(A:) = Wk, k = 1,2, . . . , n 

ierr is an error flag equal to on normal return, 

equal to i if the QR algorithm does not con- 
verge within 30 iterations on evaluating the 
ith eigenvalue, equal to -1 if n is not in range, 
and equal to -2 if one of the /3's is negative. 

The array e of dimension n is used for working space. The double precision 
routine has the name dgauss. 

We refrain here from giving numerical examples, since the use of the 
routine gauss, and the routines yet to be described, is straightforward. 

^This routine was kindly supplied to the author by Professor G.H. Golub. 
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Some use of gauss and dgauss has already been made in Examples 4.2, 4.3, 
4.4 and 5.2. 

6.2. Gauss- Radau quadrature. Wc now assume that d\ is a measure 
whose support is either bounded from below, or bounded from above, or 
both. Let, then, xq be either the infimum or the supremum of supp c?A, so 
long as it is finite. (Typically, if supp d\ = [—1, 1], then xq could be either 
-1 or +1; if supp dX = [0, oo], then xo would have to be 0; etc.). By Gauss- 
Radau quadrature we then mean a quadrature rule of maximum degree of 
exactness that contains among the nodes the point xq. It thus has the form 

n 

f{t)d\{t) = wofixo) + wkfixk) + Rnif), (6-5) 

fc=l 

and, as is well known, can be made to have degree of exactness 2n, i.e., 
Rn{f) = for all polynomials of degree < 2n. Interestingly, all nodes xq, 
xi,. . . ,Xn and weights wq, wi, . . . ,Wn can again be interpreted in terms of 
eigenvalues and eigenvectors, exactly as in the case of Gaussian quadrature 
rules, but now relative to the matrix (cf. [39]) 




VPl Oil 





g R,(«+1)X('"+1)^ 



(6.6) 

where aj- = ak{dX) {0 < k < n — 1), (3k = (3k{d^) {I < k < n) as before, but 

7r„_i(.i:o: dX) 



«n = OnidX) =Xo- Pn 



7r„(xo;o?A) 



(6.7) 



Hence, we can use the routine gauss to generate the Gauss-Radau formula. 
This is done in the following subroutine. 

subroutine radau (n , alpha , beta , end , zero , weight , ierr , e , a , b) 
dimension alpha(*) ,beta(*) ,zero(*) ,weight(*) ,e(*) ,a(*) ,b(*) 



c The arrays alpha, beta, zero, weight, e, a, b are assumed to have 
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c dimension n+1. 
c 

epsma=r Imach ( 3 ) 

c 

c epsma is the machine single precision, 
c 

npl=n+l 
do 10 k=l,npl 
a(k)=alpha(k) 
b(k)=beta(k) 
10 continue 
pO=0. 
pl=l. 

do 20 k=l,n 
pml=pO 

pO=pl 

pl= (end-a(k) ) *pO-b (k) *pml 
20 continue 

a(npl)=end-b(npl)*pO/pl 

call gauss (npl , a,b , epsma, zero , weight , ierr , e) 

return 

end 

The input variables are n, alpha, beta, end representing, respectively, n, two 
arrays of dimension n + 1 containing the ak{dX), (3k{dX), k = 0,1,2, ... ,n, 
and the "endpoint" xq. The nodes (in increasing order) of the Gauss-Radau 
formula are returned in the array zero, the corresponding weights in the 
array weight. The arrays e, a, b are working space, and ierr an error 
flag inherited from the routine gauss. The double-precision routine has the 
name dradau. 

We remark that xq could also be outside the support of dX, in which 
case the routine would generate a "Christoffel-type" quadrature rule. 

6.3. Gauss-Lobatto quadrature. Assuming now the support of dX 
bounded on either side, we let xq = inf supp (dX) and Xn+i = sup supp (dX) 
and consider a quadrature rule of the type 

„ n 

f{t)dX{t) = Wofixo) + ^ WkfiXk) + Wn+lf{Xn+l) + Rn{f) (6-8) 
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having maximum degree of exactness 2n+l. This is caUed the Gauss-Lobatto 
quadrature rule. Its nodes xq, xi, . . . , Xn+i and weights wq, wi, . . . , Wn+i 
again admit the same spectral representation as in the case of Gauss and 
Gauss-Radau formulae, only this time the matrix in question has order n+2 
and is given by (cf. [39]) 







Olr,. ^JPl 



n+1 



n+1 



a. 



n+1 



(6.9) 



Here, as before, = ak{dX) {0 < k < n), (3^ = Pk{dX) {I < k < n), and 
q:*_|_^, are the unique solution of the linear 2x2 system 

TTn+l{xo; dX) TTn{xo;dX) 
TTn+liXn+i; dX) TTniXn+r, dX) 

(6.10) 

Hence, we have the following routine for generating Gauss-Lobatto formulae: 





O^n+l 




XoTTn+lixo; dX) 








Xn+lT^n+liXn+l] dX) 



subroutine lob (n , alpha , beta , alef t , right , zero , weight , ierr , e , a , b) 
dimension alpha(*) ,beta(*) ,zero(*) ,weight(*) ,e(*) ,a(*) ,b(*) 

c 

c The arrays alpha, beta, zero, weight, e, a, b are assumed to have 

c dimension n+2. 

c 

epsma=r Imach ( 3 ) 

c 

c epsma is the machine single precision, 
c 

npl=n+l 
np2=n+2 
do 10 k=l,np2 

a(k)=alpha(k) 

b(k)=beta(k) 
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10 continue 
p01=0. 
pOr=0 . 
pll=l. 
plr=l. 

do 20 k=l,npl 
pmll=p01 
p01=pll 
piiilr=pOr 
pOr=plr 

pll= (alef t-a(k) ) *p01-b(k) *pmll 
plr= (right-a(k) ) *pOr-b(k) *pmlr 
20 continue 

det=pll*pOr-plr*p01 

a(np2) = (alef t*pll*p0r-right*plr*p01) /det 

b(np2)=(right-aleft)*pll*plr/det 

call gauss (np2 , a, b , epsma, zero , weight , ierr , e) 

return 

end 

The meaning of the input and output variables is as in the routine radau, 
the variable alef t standing for xq and right for Xn+i- The double-precision 
routine is named dlob. 

A remark analogous to the one after the routine radau applies to the 
routine lob. 
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